USE ElliDef
USE matrix_oper
USE NodeInfoDef
USE metric
USE Elliptic_Solvers
IMPLICIT NONE
TYPE (NodeInfo) :: Info


INTEGER error,N,nmax
REAL (KIND=qPrec) :: xUpper(1:3),dx,dy,dz,xmax
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl,dummy,dummy1
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:,:) :: ut,vt,wt


INTEGER i,j,k,nx,ny,nz,l
INTEGER, PARAMETER :: solver=PARDISO_SOLVER
!INTEGER, PARAMETER :: solver=IT_FOM
LOGICAL iterate









! populate Info structures

ALLOCATE(Info%ElliInfo)
xmax=3.3;nmax=5
Info%mX=(/nmax,10,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
Info%nDim=2
Info%ElliInfo%ndim=Info%ndim
Info%ElliInfo%mX=Info%mX ! last dim has to be one
!!$ALLOCATE(Info%ElliInfo%p(0:Info%mX(1)+1,0:Info%mX(2)+1,0:Info%mX(3)+1))
!!$ALLOCATE(Info%ElliInfo%pold(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)))
!!$ALLOCATE(Info%ElliInfo%qe(nx+1,ny+1,nz+1,2))
!!$ALLOCATE(Info%ElliInfo%div(nx,ny,nz))

!xUpper(1:3)=(/20.d0,30.d0,1.d0/)

!!$Info%xLower=(/-5000./2200.,0.,0./)
!!$
!!$Info%xUpper=(/15000./2200,1.,1./)

Info%xLower=(/-xmax,0.d0,0.d0/)
!!$
Info%xUpper=(/xmax,1.d0,1.d0/)


Info%dX(1:2)=(Info%xUpper(1:2)-Info%xLower(1:2))/Info%mX(1:2);Info%dX(3)=1
Info%ElliInfo%dX(1:2)=Info%dX(1:2)
Info%ElliInfo%dX(3)=Info%dX(3) ! just to have it set to something
Info%dmax=0.

!Info%mthbc=InternalBC ! assume all Neumann conditions
Info%ElliInfo%mthbc(1:6)=0
Info%ElliInfo%mthbc(4)=0*InternalBC




CALL setprob(Info)
PRINT*,'Dimension of problem',Info%mX

ALLOCATE(ut(1:ny,1:nz,2),vt(1:nx,1:nz,2),wt(nx,ny,2))
CALL DEFINE_TOPO(Info) ! define the quantities used to calculate the metric
error=CREATE_OPERATORS_THIN(Info%ElliInfo)

! zero BC here to be safe 
DO i=1,2*Info%nDim
   Info%ElliInfo%bc(i)%p=4*one
END DO

IF(.NOT.check(Info%ElliInfo,CHECK_BC)) THEN
   PRINT*,'BC are not consistent'
   STOP
END IF






! calculates divergence


u=>Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);


CALL RANDOM_NUMBER(u)
CAlL RANDOM_NUMBER(v)
!!$u=zero;v=zero
!!$v(:,:,1)=SPREAD(-one*.5*(cos(8.*atan(1.d0)*((/1:nx/)-0.5)/float(nx))),DIM=2,NCOPIES=ny+1)
ut=zero;vt=zero
CALL RANDOM_NUMBER(ut)
CALL RANDOM_NUMBER(vt)
ut(:,1,1)=ut(:,1,1)-SUM(ut(:,1,1))/ny
ut(:,1,2)=ut(:,1,2)-SUM(ut(:,1,2))/ny
vt(:,1,1)=vt(:,1,1)-SUM(vt(:,1,1))/nx
vt(:,1,2)=vt(:,1,2)-SUM(vt(:,1,2))/nx


CALL take_div(Info%Elliinfo,.false.)
PRINT*,Info%ElliInfo%qe(3,:,1,2)
PRINT*,'**************'
        PRINT*,SUM(Info%ElliInfo%div**2)



! only need OP= optional to use IT wth PARDISO precon
i=POISSON_SOLVERS_THIN(Info%ElliInfo,INIT)!,OP=CSR_PRECOND)

CALL hodge_thin(Info%ElliInfo,300,.false.,ut,vt)


CALL take_div(Info%Elliinfo,.false.)
PRINT*,'****************'

PRINT*,Info%ElliInfo%p(3,:,1)-4.

i=POISSON_SOLVERS_THIN(Info%ElliInfo,FINISH)


PRINT*,'parison_mopup'


10 open(10,file='u.dat',form='unformatted')
write(10) Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1)
close(10)
open(10,file='w.dat',form='unformatted')
write(10) Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2)
close(10)
!!$open(10,file='p.dat',form='unformatted')
!!$write(10) wt(0:nx+1,0:ny+1)
!!$close(10)
open(10,file='div.dat',form='unformatted')
write(10) Info%ElliInfo%div(1:nx,1:ny,1:nz)
close(10)
PRINT*,SUM(Info%ElliInfo%div**2)
error=FREEUP_OPERATORS(Info%ElliInfo)
 
DEALLOCATE(Info%ElliInfo)
DEALLOCATE(Info%topo%h,Info%topo%hx)
DEALLOCATE(Info%topo)
DEALLOCATE(ut,vt)


END

